Making a map with points in R with leaflet

We're going to be doing a lot of the stuff from the RStudio leaflet tutorial. This walkthrough is a variation of Amelia McNamara's version.

Load leaflet first:

library(leaflet)

Let's start with a simple map

m <- leaflet() %>%
  setView(-72.518978, 42.381050, zoom = 1) %>%
  addTiles() %>%  # Load the default OpenStreetMap tile set (images that make up the map)
  addMarkers(lng=-72.518978, lat=42.381050, popup="DSDP!")

m

Note that the syntax here is similar to that of ggplot2. Please consult the Leaflet documentation for more details.

Okay, let's do something with actual data

Lets look at storm data from the NOAA. It comes in a few files that we need to join together in order to use. For convenience, I've put them on the workshop website. You can download them manually, or you can use the getURL() function from RCurl to grab them from the web as I do below (they're big, so it takes a sec):

library(readr)
library(dplyr)
library(RCurl)

# One file contains information about the storm
stormdetails <- getURL("https://jcrouser.github.io/MassMutual-DataVis/datasets/StormEvents_details-ftp_v1.0_d2016_c20160810.csv")
stormdetails <- read_csv(stormdetails)
## Warning: 176 parsing failures.
##  row                col           expected actual         file
## 1233 TOR_OTHER_WFO      1/0/T/F/TRUE/FALSE   TAE  literal data
## 1233 TOR_OTHER_CZ_STATE 1/0/T/F/TRUE/FALSE   AL   literal data
## 1233 TOR_OTHER_CZ_FIPS  1/0/T/F/TRUE/FALSE   045  literal data
## 1233 TOR_OTHER_CZ_NAME  1/0/T/F/TRUE/FALSE   DALE literal data
## 1498 TOR_OTHER_WFO      1/0/T/F/TRUE/FALSE   TAE  literal data
## .... .................. .................. ...... ............
## See problems(...) for more details.
# The other contains actual location data
stormlocs <- getURL("https://jcrouser.github.io/MassMutual-DataVis/datasets/StormEvents_locations-ftp_v1.0_d2016_c20160810.csv")
stormlocs <- read_csv(stormlocs)

# We'll want to use them together, so we'll use a join
storms <- stormlocs %>%
  left_join(stormdetails, by="EVENT_ID")
glimpse(storms)
## Rows: 10,818
## Columns: 61
## $ YEARMONTH          <dbl> 201603, 201603, 201603, 201603, 201603, 201603, 201603, 201603, 201603, 2…
## $ EPISODE_ID.x       <dbl> 102667, 102667, 102667, 102667, 102693, 102693, 102693, 102693, 102738, 1…
## $ EVENT_ID           <dbl> 615187, 615188, 615189, 615190, 625566, 625566, 625567, 625567, 613936, 6…
## $ LOCATION_INDEX     <dbl> 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3…
## $ RANGE              <dbl> 2.82, 1.84, 0.00, 0.00, 3.00, 5.00, 0.00, 0.50, 0.91, 0.86, 5.03, 5.00, 5…
## $ AZIMUTH            <chr> "SSE", "WSW", "N", "N", "ESE", "SW", "N", "E", "SSW", "NW", "S", "S", "S"…
## $ LOCATION           <chr> "ANTIOCH", "SCULLIN", "CONNERVILLE", "ATOKA", "GASTON", "BEAVERTON", "NEH…
## $ LATITUDE           <dbl> 34.6800, 34.5100, 34.4500, 34.3800, 45.4134, 45.4288, 45.7200, 45.7200, 3…
## $ LONGITUDE          <dbl> -97.4000, -96.8900, -96.6300, -96.1300, -123.0729, -122.8930, -123.9000, …
## $ LAT2               <dbl> 3440800, 3430600, 3427000, 3422800, 4524804, 4525728, 4543200, 4543200, 3…
## $ LON2               <dbl> 9724000, 9653400, 9637800, 967800, 1234374, 12253580, 12354000, 12353376,…
## $ BEGIN_YEARMONTH    <dbl> 201603, 201603, 201603, 201603, 201603, 201603, 201603, 201603, 201603, 2…
## $ BEGIN_DAY          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7…
## $ BEGIN_TIME         <dbl> 2, 52, 110, 155, 1630, 1630, 1730, 1730, 1621, 135, 520, 520, 520, 520, 5…
## $ END_YEARMONTH      <dbl> 201603, 201603, 201603, 201603, 201603, 201603, 201603, 201603, 201603, 2…
## $ END_DAY            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7…
## $ END_TIME           <dbl> 2, 52, 110, 155, 1645, 1645, 1745, 1745, 1621, 135, 526, 526, 526, 526, 5…
## $ EPISODE_ID.y       <dbl> 102667, 102667, 102667, 102667, 102693, 102693, 102693, 102693, 102738, 1…
## $ STATE              <chr> "OKLAHOMA", "OKLAHOMA", "OKLAHOMA", "OKLAHOMA", "OREGON", "OREGON", "OREG…
## $ STATE_FIPS         <dbl> 40, 40, 40, 40, 41, 41, 41, 41, 22, 5, 53, 53, 53, 53, 53, 53, 53, 53, 53…
## $ YEAR               <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2…
## $ MONTH_NAME         <chr> "March", "March", "March", "March", "March", "March", "March", "March", "…
## $ EVENT_TYPE         <chr> "Hail", "Hail", "Hail", "Hail", "Thunderstorm Wind", "Thunderstorm Wind",…
## $ CZ_TYPE            <chr> "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C"…
## $ CZ_FIPS            <dbl> 49, 99, 69, 5, 67, 67, 57, 57, 59, 99, 47, 47, 47, 47, 19, 19, 19, 19, 65…
## $ CZ_NAME            <chr> "GARVIN", "MURRAY", "JOHNSTON", "ATOKA", "WASHINGTON", "WASHINGTON", "TIL…
## $ WFO                <chr> "OUN", "OUN", "OUN", "OUN", "PQR", "PQR", "PQR", "PQR", "SHV", "SHV", "OT…
## $ BEGIN_DATE_TIME    <chr> "01-MAR-16 00:02:00", "01-MAR-16 00:52:00", "01-MAR-16 01:10:00", "01-MAR…
## $ CZ_TIMEZONE        <chr> "CST-6", "CST-6", "CST-6", "CST-6", "PST-8", "PST-8", "PST-8", "PST-8", "…
## $ END_DATE_TIME      <chr> "01-MAR-16 00:02:00", "01-MAR-16 00:52:00", "01-MAR-16 01:10:00", "01-MAR…
## $ INJURIES_DIRECT    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ INJURIES_INDIRECT  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ DEATHS_DIRECT      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ DEATHS_INDIRECT    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ DAMAGE_PROPERTY    <chr> "0.00K", "0.00K", NA, "0.00K", "3.00K", "3.00K", "24.00K", "24.00K", "0.0…
## $ DAMAGE_CROPS       <chr> "0.00K", "0.00K", NA, "0.00K", "0.00K", "0.00K", "0.00K", "0.00K", "0.00K…
## $ SOURCE             <chr> "Trained Spotter", "Amateur Radio", "Emergency Manager", "Public", "Publi…
## $ MAGNITUDE          <dbl> 1.00, 1.00, 1.75, 1.00, 47.00, 47.00, 47.00, 47.00, 1.00, 52.00, NA, NA, …
## $ MAGNITUDE_TYPE     <chr> NA, NA, NA, NA, "EG", "EG", "EG", "EG", NA, "EG", NA, NA, NA, NA, NA, NA,…
## $ FLOOD_CAUSE        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Heavy Rain", "Heavy Rain", "Heav…
## $ CATEGORY           <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_F_SCALE        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_LENGTH         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_WIDTH          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_OTHER_WFO      <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_OTHER_CZ_STATE <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_OTHER_CZ_FIPS  <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ TOR_OTHER_CZ_NAME  <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ BEGIN_RANGE        <dbl> 3, 2, 0, 0, 3, 3, 0, 0, 1, 1, 5, 5, 5, 5, 10, 10, 10, 10, 1, 1, 1, 1, 1, …
## $ BEGIN_AZIMUTH      <chr> "SSE", "WSW", "N", "N", "ESE", "ESE", "N", "N", "SSW", "NW", "S", "S", "S…
## $ BEGIN_LOCATION     <chr> "ANTIOCH", "SCULLIN", "CONNERVILLE", "ATOKA", "GASTON", "GASTON", "NEHALE…
## $ END_RANGE          <dbl> 3, 2, 0, 0, 5, 5, 1, 1, 1, 1, 5, 5, 5, 5, 10, 10, 10, 10, 1, 1, 1, 1, 1, …
## $ END_AZIMUTH        <chr> "SSE", "WSW", "N", "N", "SW", "SW", "E", "E", "SSW", "NW", "S", "S", "S",…
## $ END_LOCATION       <chr> "ANTIOCH", "SCULLIN", "CONNERVILLE", "ATOKA", "BEAVERTON", "BEAVERTON", "…
## $ BEGIN_LAT          <dbl> 34.6800, 34.5100, 34.4500, 34.3800, 45.4134, 45.4134, 45.7200, 45.7200, 3…
## $ BEGIN_LON          <dbl> -97.4000, -96.8900, -96.6300, -96.1300, -123.0729, -123.0729, -123.9000, …
## $ END_LAT            <dbl> 34.6800, 34.5100, 34.4500, 34.3800, 45.4288, 45.4288, 45.7200, 45.7200, 3…
## $ END_LON            <dbl> -97.4000, -96.8900, -96.6300, -96.1300, -122.8930, -122.8930, -123.8896, …
## $ EPISODE_NARRATIVE  <chr> "A line of storms formed along a boundary in the Texas panhandle on the e…
## $ EVENT_NARRATIVE    <chr> "Hail covered the entire ground with some hail drifts.", NA, "No damage w…
## $ DATA_SOURCE        <chr> "CSV", "CSV", "CSV", "CSV", "CSV", "CSV", "CSV", "CSV", "CSV", "CSV", "CS…

Don't worry if you don't really understand the syntax above. Next week, Ben Baumer will teach you all about SQL and joins.

Mapping lightning strikes

Let's pull out the lightning strikes and map them as part of a dplyr chain

lightning_map <- storms %>%               # Start with the storm data
  filter(EVENT_TYPE == "Lightning") %>%   # Filter down to just the lightning events
  leaflet() %>%                           # Pipe into a leaflet map
    addMarkers(~LONGITUDE, ~LATITUDE, popup = "Zap!") %>%   # Add markers at each LONGITUDE/LATITUDE pair
    addProviderTiles("Stamen.Toner")        # We'll use black/white tiles for dramatic effect

lightning_map

Challenge:

  • Find another storm type to map
  • Bonus-- add popups!

One approach

mtw_map <- storms %>%
  filter(EVENT_TYPE == "Marine Thunderstorm Wind") %>%
  leaflet() %>%
    addProviderTiles("Stamen.Toner") %>% 
    addMarkers(~LONGITUDE, ~LATITUDE, popup = ~EVENT_NARRATIVE)

mtw_map

Polygons

First let's grab some data (and do a little conversion)

tornados <- storms %>%
  filter(EVENT_TYPE=="Tornado") %>%
  mutate(DAMAGE_PROPERTY = as.numeric(sub("K", "", DAMAGE_PROPERTY, fixed = TRUE)))
## Warning in mask$eval_all_mutate(dots[[i]]): NAs introduced by coercion

Then let's start with something easy(ish) -- circles

m <- leaflet(data = tornados) %>%
  addProviderTiles("Stamen.Toner") %>% 
  addCircles(~LONGITUDE, ~LATITUDE, 
             weight = 1, 
             radius = ~DAMAGE_PROPERTY*100, # Map the radius of the circle to amount of damage
             popup = ~EVENT_NARRATIVE)      # Include details about the tornado
m

Polgyons come in shapefiles

Most boundaries (state, national, etc) are provided in terms of polygons. Major mapping software ArcGIS, from ESRI, has essentially set the standard formats. There are many files with different extensions: .prj (the projection), .shp (the shapefile), .cpg (??), .dbf (??), .shx (??).

You need special software or packages to work with shapefiles.

State shapefiles

I got these from the Census. You can choose the resolution.

If you want, the zipfile of the shapes I used is here.

We're going to use the maptools package to deal with shapefiles.

library(maptools)

# Remember when we talked about map projections?
crswgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Now we can load the shapefile using the correct projection
states = readShapePoly("datasets/cb_2015_us_state_500k/cb_2015_us_state_500k.shp",
                       proj4string = crswgs84,
                       verbose = TRUE)
## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
## Shapefile type: Polygon, (5), # of Shapes: 56

In RStudio, you can click on the states object to see what a Large SpatialPolygonsDataFrame looks like. It contains both @data and @polygons.

Exploring shapefiles

Let's start by mapping some boring, internal data from shapefile about the amount of water per state

states %>%
  leaflet() %>%
  setView(-95.976807, 40.829587, zoom = 3) %>%
  addProviderTiles("Stamen.Toner") %>%
  addPolygons(stroke = FALSE, 
              fillOpacity = 0.5, 
              smoothFactor = 0.5, 
              color = ~colorQuantile("BrBG", states$AWATER)(AWATER)
  )

Challenge

Put it all together! See if you can make a map of your own data, and embed it within a dashboard. Bonus: can you make them interact with one another?